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Xfy Abstract 

0^ We show how one may analytically compute the stationary density of the distribution of molecular constituents in populations of 
cells in the presence of noise arising from either bursting transcription or translation, or noise in degradation rates arising from low 

! ! numbers of molecules. We have compared our results with an analysis of the same model systems (either inducible or repressible 

^ | operons) in the absence of any stochastic effects, and shown the correspondence between behaviour in the deterministic system and 
m 'the stochastic analogs. We have identified key dimensionless parameters that control the appearance of one or two steady states in 
<— 5 the deterministic case, or unimodal and bimodal densities in the stochastic systems, and detailed the analytic requirements for the 
Q occurrence of different behaviours. This approach provides, in some situations, an alternative to computationally intensive stochas- 
tic simulations. Our results indicate that, within the context of the simple models we have examined, bursting and degradation noise 
cannot be distinguished analytically when present alone. 
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1. Introduction 

In neurobiology, when it became clear that some of the fluc- 
tuations seen in whole nerve recording, and later in single cell 
recordings, were not simply measurement noise but actual fluc- 
tuations in the system being studied, researchers very quickly 
started wondering to what extent these fluctuations actually 
played a role in the operation of the nervous system. 
, Much the same pattern of development has occurred in cel- 
lular and molecular biology as experimental techniques have 
allowed investigators to probe temporal behaviour at ever finer 
levels, even to the level of individual molecules. Experimental- 
ists and theoreticians alike who are interested in the regulation 
of gene networks are increasingly focussed on trying to access 
the role of various types of fluctuations on the operation and 
fidelity of both simple and complex gene regulatory systems. 
Recent reviews dKaernet all 120051; [Raj and van Qudenaardeni 
l2008t IShah rezaei and SwainL l2008bl " give an interesting per- 
spective on some of the issues confronting both experimental- 
ists and modelers. 



Typically, the discussion seems to focus on whether fluc- 
tuations can be considered as extrinsic to the system under 



consideration (Shahrezaei et al., 2008; Ochab-Marcinek, 2008 



201 Oh . or whether they are an intrinsic part of the fundamen- 
tal processes they are affecting (e.g. bursting , see below) 



Elowitz et al 



The d ichotomy is rarely so sharp however, but 
(120021) have used an elegant expe rimental technique to dis 
tingui sh between t he two , se e also Raser and O'Sheal d2004l) . 



while ISwain et all J2002h and lScott et all d2006h have laid the 
groundwork for a theoretical consideration of this question. 
One issue that is raised persistently in considerations of the 
role of fluctuations or noise in the operation of ge ne regulator 
netwo rks is whether or not they are "beneficial" dBlake et 



ory 
al. 
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20061) or "detrimental" dFraser et all 120041) to the operation of 
the system under consideration. This is, of course, a question 
of definition and not one that we will be further concerned with 
here. 



Here, we consider in detail the density of the molecular 
distributions in generic bacterial operons in the presence of 
'bursting' (commonly known as intrinsic noise in the bio- 
logical literature) as well as inherent (extrinsic) noise using 
an analytical approach. Our work is motivated by the well 
documented production of mRNA and/or protei n in stochas 
tic bursts in both prokaryotes and eukaryotes ( Blake et al 



20031; ICai et all 120061: IChubb et all 120061: iGolding et all 12005 
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Rai et all 120061: ISieal et all 120061: lYu et all 120061). and fol- 
lows o t her contribution s by, f o r example, Kep ler and Elston 



(l200ll), iFriedman et al l d2006k . iBobrowski et alj d2007l) and 



Shahrezaei and Swainl d2008al) . 

In Section [2] we develop the concept of the operon and treat 
simple models of the classic inducible and repressible operon. 
Section|4]considers the effects of bursting alone in an ensemble 
of single cells. Section [5] then examines the situation in which 
there are continuous white noise fluctuations in the dominant 
species degradation rate in the absence of bursting. 



permease, and thiogalactoside transacetylase respectively. The 
enzyme /3-galactosidase is active in the conversion of lactose 
into allolactose and then the conversion of allolactose into glu- 
cose. The lac permease is a membrane protein responsible for 
the transport of extracellular lactose to the interior of the cell. 
(Only the transacetylase plays no apparent role in the regula- 
tion of this system.) The regulatory gene lad, which is part of 
a different operon, codes for the lac repressor, which is trans- 
formed to an inactive form when bound with allolactose, so in 
this system allolactose functions as the effector molecule. 



2. Generic operons 

2.1. The operon concept 

The so-called 'central dogma' of molecular biology is sim- 
ple to state in principle, but complicated in its detail. Namely 
through the process of transcription of DNA, messenger RNA 
(mRNA, M) is produced and, in turn, through the process of 
translation of the mRNA proteins (intermediates, /) are pro- 
duced. There is often feedback in the sense that molecules (en- 
zymes, E) whose production is controlled by these proteins can 
modulate the translation and/or transcription processes. In what 
follows we will refer to these molecules as effectors. We now 
consider both the transcription and translation process in more 
detail. 

In the transcription process an amino acid sequence in the 
DNA is copied by the enzyme RNA polymerase (RNAP) to 
produce a complementary copy of the DNA segment encoded 
in the resulting RNA. Thus this is the first step in the transfer 
of the information encoded in the DNA. The process by which 
this occurs is as follows. 

When the DNA is in a double stranded configuration, the 
RNAP is able to recognize and bind to the promoter region of 
the DNA. (The RNAP/double stranded DNA complex is known 
as the closed complex.) Through the action of the RNAP, the 
DNA is unwound in the vicinity of the RNAP/DNA promoter 
site, and becomes single stranded. (The RNAP/single stranded 
DNA is called the open complex.) Once in the single stranded 
configuration, the transcription of the DNA into mRNA com- 
mences. 

In prokaryotes, translation of the newly formed mRNA com- 
mences with the binding of a ribosome to the mRNA. The 
function of the ribosome is to 'read' the mRNA in triplets of 
nucleotide sequences (codons). Then through a complex se- 
quence of events, initiation and elongation factors bring transfer 
RNA (tRNA) into contact with the ribosome-mRNA complex to 
match the codon in the mRNA to the anti-codon in the tRNA. 
The elongating peptide chain consists of these linked amino 
acids, and it starts folding into its final conformation. This fold- 
ing continues until the process is complete and the polypeptide 
chain that results is the mature protein. 

The lactose (lac) operon in bacteria is the paradigmatic ex- 
ample of this concept and this much studied system consists of 
three structural genes named lacZ, lacY, and lack. These three 
genes contain the code for the ultimate production, through the 
translation of mRNA, of the intermediates yS-galactosidase, lac 



2.2. The transcription rate function 

In this section we examine the molecular dynamics of both 
the classical inducible and repressible operon to derive expres- 
sions for the dependence of the transcription rate on effector 
levels. (When the transcription rate is constant and indepen- 
dent of the effector levels we will refer to this as the no control 
situation.) 

2.2.1. Inducible regulation 

For a typical inducible regulatory situation (such as the lac 
operon), in the presence of the effector molecule the repressor 
is inactive (is unable to bind to the operator region preceding 
the structural genes), and thus DNA transcription can proceed. 
Let R denote the repressor, E the effector molecule, and O the 
operator. The effector is known to bind with the active form R 
of the repressor. We assume that this reaction is of the form 



R + nE% RE„ 



Ki = 



RE„ 
RE"' 



(1) 



where n is the effective number of molecules of effector re- 
quired to inactivate the repressor R. Furthermore, the operator 
O and repressor R are assumed to interact according to 



O + R^OR 



K 2 



OR 
OR' 



Let the total operator be O tol : 

0,„, = + OR = + K 2 R = 0(l+ K 2 R), 

and the total level of repressor be R tot : 

R, ot = R + K l RE" +K 2 OR. 

The fraction of operators not bound by repressor (and therefore 
free to synthesize mRNA) is given by 



f(E) 



O 



1 



0, , 1 + K 2 R 

If the amount of repressor R bound to the operator O is small 
R tot ^R + KiR ■ E" = R(l + K\E") 



R 



R„ 



1+KiE" 



2 



and consequently 

/(£) = 



1 +KiE" _ 1 +K 1 E n 
1 +K 2 R,ot +K l E" ~ K + K\E" ' 



(2) 



where K — 1 + K 2 R, ot . There will be maximal repression when 
E — but even then there will still be a basal level of mRNA 
production proportional to K~ l (which we call the fractional 
leakage). 

If the maximal DNA transcription rate is tp m (in units of in- 
verse time) then, under the assumption that the rate of transcrip- 
tion <p in the entire population is proportional to the fraction / 
of unbound operators, the variation <p of the DNA transcription 
rate with the effector level is given by <p = <p m f, or 



<p(E) = tp„ 



1 + KiE" 
K + KiE" 



(3) 



2.2.2. Repressible regulation 

In the classic example of a repressible system (such as the trp 
operon) in the presence of the effector molecule the repressor is 
active (able to bind to the operator region), and thus block DNA 
transcription. We use the same notation as before, but now note 
that the effector binds with the inactive form R of the repressor 
so it becomes active. We assume that this reaction is of the same 
form as in Equation [TJ The difference now is that the operator 
O and repressor R are assumed to interact according to 

„ ORE,, 
O+R E" ^ ORE,, K 2 = 



ORE,, 
The total operator is now given by 

O tot = + ORE,, = O + K 2 R E" = 0(1 + K 2 R ■ E"), 
so the fraction of operators not bound by repressor is given by 

1 



f(E) 



O 
~0~o< 



1 + K-,R ■ E" 



Again assuming that the amount of repressor R bound to the 
operator O is small we have 



f(E) 



1 +KiE" 



1 +K\E n 



l+(Ki+ K 2 R I0I )E" 1 + KE" 



where K — K\ + K 2 R„„. Now there will be maximal repression 
when E is large, but even at maximal repression there will still 
be a basal level of mRNA production proportional to K\K X < 
1. The variation of the DNA transcription rate with effector 
level is given by (p — <p m f or 



<p(E) = (p„ 



1 +KiE" 



1 + KE" 

Both (0) and (0| are special cases of the function 
1 + K X E" 



* (E) = ^A + BE» 
where A, B > are given in Table [TJ 



= tPmf(E). 



(4) 



(5) 



parameter 


inducible 


repressible 


A 


K = 1 + K 2 R lo , 


1 


B 


Ki 


K = K\ + K 2 R, i 


B 

A 


Ki 
K 


A 


A — A 
A — /i 


A 


1 
1 


A = bk- 1 


l 


KK- 1 


nA\ A) 


*«. K - l >o 

n K 


*.*-* <0 
n K 



Table 1 : Definitions of the parameters A, B, A, A and 8. See the text and Section 
I2.2l for more detail. 



2.3. Deterministic operon dynamics in a population of cells 



The reader may wish to consult lPolvnikis et al. (2009) for an 
interesting survey of techniques applicable to this approach. 

We first consider a large population of cells, each of which 
contains one copy of a particular operon, and let (M, /, E) de- 
note mRNA, intermediate protein, and effector levels respec- 
tively in the population. Then for a generic operon with a max- 
imal level of transcription (in co ncentration unit s), we have 
dynamics described b y the system ((Griffith, ll968ajbL Othmer, 
1 1976tlSelgradel.il 9791) 



dM - 

— = b d tp„,f(E) - y M M, 
dt 

dl 

dE 

-=B E I-y E E. 



(6) 
(7) 
(8) 



Here we assume that the rate of mRNA production is propor- 
tional to the fraction of time the operator region is active, and 
that the rates of intermediate and enzyme production are simply 
proportional to the amount of mRNA and intermediate respec- 
tively. All three of the components (M, I, E) are subject to ran- 
dom loss. The function / is calculated in the previous section. 

It will greatly simplify matters to rewrite Equations |6][8] by 
defining dimensionless concentrations. To this end we first 
rewrite Equation[5]in the form 

<p(e) = ip m f(e), (9) 

where <p m (dimensionless) is defined by 

<f m B E Bi l+e" 

<Pm = and /(e) = — — — , (10) 

JmJeJi A + Ae" 

A and A are defined in Table [TJ and we have defined a dimen- 
sionless effector concentration (e) through 



E -r\e with r\ 



1 



Further defining dimensionless intermediate (/) and mRNA 
concentrations (m) through 



in — and M = mil , 

'p E 'PePi 



Equations [6]|8]can be written in the equivalent form 



dm 

— = jM[Kdf{e)-m\, 
at 



di 
dt 

de 

~di 



ji(m - i), 
7eU ~ e), 



where 



Kd = b d <p m and b d 



(ID 



are dimensionless constants. 

For notational simplicity, henceforth we denote dimension- 
less concentrations by (m,i,e) = (xi, Xz, X3), and subscripts 
(M, I,E) = (1,2, 3). Thus we have 



-JT = rdKdfixi) - xi], 
dt 



dX2 

It 



dx 3 

—77 = 73(X2 - X3). 

dt 



(12) 
(13) 
(14) 



In each equation, y, for i = 1, 2, 3 denotes a net loss rate (units 
of inverse time), and thus Equations [T2lfT4l are not in dimen- 
sionless form. 

The dynamics of this classic operon model can be fully ana- 
lyzed. Let X = (jci,jc2i x i) and denote by S,(X) the flow gener- 
ated by the system (TT2T>-(fT4b. For both inducible and repressible 
operons, for all initial conditions X° = (x^, x^, x?) e the flow 
S,(X°) eM+forf > 0. 

Steady states of the system ([T2li-(fT4li are in a one to one cor- 
respondence with solutions of the equation 



-=/(*) 

Kd 



(15) 



and for each solution x* of Equation Q3] there is a steady state 
X* = (x\,x* v x*) of (fl2ll-(fT4l given by 

Whether there is a single steady state X* or there are multiple 
steady states will depend on whether we are considering a re- 
pressible or inducible operon. 

2.3.1. No control 

In this case, f(x) = 1, and there is a single steady state x* = 
Kd that is globally asymptotically stable. 




Figure 1 : Schematic illustration of the possibility of one, two or three solu- 
tions of Equation 1 1 5 1 for varying values of kj with inducible regulation. The 
monotone increasing graph is the function / of Equation 1101 and the straight 
lines correspond to x/kj for (in a clockwise direction) kj £ [0, kj ), = kj-, 
Kd £ (Kd-, kj = KJ+, and kj + < kj. This figure was constructed with n = 4 
and K = 10 for which kj- = 3.01 and kj + = 5.91 as computed from 4 1 8i . See 
the text for further details. 



2.3.2. Inducible regulation 

Single versus multiple steady states. For an inducible operon 
with / given by Equation [2] there may be one (X^ or X^), two 
(X* V X* = X* or X* = X*,X*), or three (X*,X*,X*) steady 
states, with the ordering < X* < X* 2 < XX, corresponding 
to the possible solutions of Equation [T31 (cf. Figure [TJ. The 
smaller steady state (X*) is typically referred to as an uninduced 
state, while the largest steady state (X!) is called the induced 
state. The steady state values of x are easily obtained from ( IT5l ) 
for given parameter values, and the dependence on K d for n — 4 
and a variety of values of K is shown in Figure Q] Figure [2] 
shows a graph of the steady states x* versus K d for various val- 
ues of the leakage parameter K. 

Analytic conditions for the existence of one or more steady 
states can be obtained by using Equation[l5]in conjunction with 
the observation that the delineation points are marked by the 
values of Kd at which x/Kd is tangent to f(x) (see Figure [TJ. 
Simple differentiation of <TT~5T > yields the second condition 



1 



(16) 



K d n{K-V) (K + x") 2 ' 

From equations (fl3T l and ( fTol l we obtain the values of x at which 
tangency will occur: 



x± = 





K + l 






n 






K-l 





(17) 



The two corresponding values of K d at which a tangency occurs 
are given by 

K + x'i 

Kd± = **-. — r- ( 18 > 

1 + x% 

(Note the deliberate use of x T as opposed to x ± .) 



4 



K 



Figure 2: Full logarithmic plot of the steady state values of x* versus 
for an inducible system, obtained from Equation [15] for n = 4 and K = 
2,5, 10, and 15 (left to right) illustrating the dependence of the occurrence of 
bistability on K. See the text for details. 



A necessary condition for the existence of two or more steady 
states is obtained by requiring that the square root in $TT\ be 
non-negative, or 

(19) 



K > 



From this a second necessary condition follows, namely 



n + 1 J n + 1 



1 



1 



(20) 



Further, from Equations[l5]and[T6]we can delineate the bound- 
aries in {K,Kd) space in which there are one or three locally 
stable steady states as shown in Figure[3] There, we have given 
a parametric plot (x is the parameter) of k c \ versus K, using 



K(x) 



x"[x" + (n + l)] 
(n - \)x n - 1 



and Kd(x) 



[K(x) + x"] 2 
nx"- l [K(x)- 1]' 



for n — 4 obtained from Equations Q3] and Q~6] As is clear from 
the figure, when leakage is appreciable (small K, e.g for n = 4, 
K < (5/3) 2 ) then the possibility of bistable behaviour is lost. 

Remark 1. Some general observations on the influence of n, 
K, and kj on the appearance of bistability in the deterministic 
case are in order. 

1 . The degree of cooperativity (n) in the binding of effector 
to the repressor plays a significant role. Indeed, n > 1 is a 
necessary condition for bistability. 

2. If n > 1 then a second necessary condition for bistability 
is that K satisfies Equation [TP] so the fractional leakage 
(K^ 1 ) is sufficiently small. 

3. Furthermore, Kd must satisfy Equation\20\which is quite 
instructive. Namely for n — > oo the limiting lower limit is 
K,j > 1 while for n — > 1 the minimal value of Kd becomes 
quite large. This simply tells us that the ratio of the product 



Figure 3: In this figure we present a parametric plot (for n = 4) of the bifur- 
cation diagram in (K,kj) parameter space delineating one from three steady 
states in a deterministic inducible operon as obtained from Equations 1151 and 
1161 The upper (lower) branch corresponds to Kd- (#</+), and for all values 
of (K,K,j) in the interior of the cone there are two locally stable steady states 
X*, Xj, while outside there is only one. The tip of the cone occurs at (K, kj) = 
((5/3) 2 , (5/3)^573) as given by Equations[l9]and[20] For K e [0, (5/3) 2 ) there 
is but a single steady state. 



of the production rates to the product of the degradation 
rates must always be greater than 1 for bistability to occur, 
and the lower the degree of cooperativity («) the larger the 
ratio must be. 

4. Ifn, K and Kd satisfy these necessary conditions then bista- 
bility is only possible if Kd 6 [Kd-,Kd+] (c.f Figure\3}. 

5. The locations of the minimal (x_) and maximal (x + ) values 
of x bounding the bistable region are independent of Kd- 

6. Finally 

(a) (x + — X-) is a decreasing function of increasing nfor 
constant Kd, K 

(b) (x + — X-) is an increasing function of increasing K 
for constant n, Kd. 

Local and global stability. The local stability of a steady state 
x* is determined by the solutions of the eigenvalue equation 
dYildirim et al.U2004 



(A + y 1 )(A + y 2 )(A + y 3 )- 7m y3Kdf: = 0, ft = f'(x). (21) 
Set 



3 3 3 

= ^ Ju a z = ^ jijj, a 3 = (1 - K d fl) Y\ 7u 



1=1 



so (I2TI 1 can be written as 



A 3 + a\A 2 + a z A + 03 = 0. 



(22) 



By Descartes's rule of signs, (l22l will have either no positive 
roots for /„' e [0, k^) or one positive root otherwise. With 
this information and using the notation SN to denote a locally 
stable node, HS a half or neutrally stable steady state, and US 
an unstable steady state (saddle point), then there will be: 



5 




0.2 0.4 0.6 0.8 



Figure 4: Schematic illustration that there is only a single solution of Equation 
1151 for all values of kj with repressible regulation. The monotone decreasing 
graph is / for a repressible operon, while the straight lines are x/kj . This figure 
was constructed with n = 4 and A = 10. See the text for further details. 



• A single steady state X* { (SN), for kj e [0, Kd-) 

• Two coexisting steady states X* (SN) and Xi = X^ (HS, 
born through a saddle node bifurcation) for K d = K d - 

• Three coexisting steady states X\{SN),X* 2 {US),Xl (SN) 
for K d e (K d -,K d+ ) 

• Two coexisting steady states X* = X* 2 (HS at a saddle node 
bifurcation), and X^ (SN) for K d = K d + 

• One steady state X^ (SN) for K d + < Kd. 

For the inducible operon, other work extends these local sta- 
bility considerations and we have the following result charac- 
terizing the global behaviour: 



Theorem 1. AOthmer \Smitn 179951 Proposition 2.1, 

Chapter 4) For an inducible operon with tp given by Equation\3\ 
define I\ — [l/K, 1]. There is an attracting box Bj C defined 
by 

Bi — {(x\, xi, x-y) : x; e //, i — 1,2,3} 

such that the flow S , is directed inward everywhere on the sur- 
face ofBj. Furthermore, allX* e Bi and 

1. If there is a single steady state, i.e. X* { for Kd 6 [0, Kd-), or 
X^for Kd+ < Kd, then it is globally stable. 

2. If there are two locally stable nodes, i.e. X* and XZ for 
Kd £ (Kd-, K d+), then all flow s S (X°) are attracted to one of 
them. (See Selsrade ( 197$il for a delineation of the basin 
of attraction ofX* and Xy) 

2.3.3. Repressible regulation 

As illustrated in Figure |4] the repressible operon has a single 
steady state corresponding to the unique solution x* of Equation 
[T31 To determine its local stability we apply the Routh-Hurwitz 



criterion to the eigenvalue equation (12211 . The steady state cor- 
responding to x* will be locally stable (i.e. have eigenvalues 
with negative real parts) if and only if a\ > (always the case) 
and 

aifl2 - «3 > 0. (23) 

The well known relation between the arithmetic and geometric 
means 



1 " 
n ■f-H 



v;=i 



l/n 



when applied to both a i and «2 gives, in conjunction with Equa- 
tion [23] 

3 

a v a 2 - «3 > (8 + K d fl) Y\ Ji > 0. 

Thus as long as / s ' > -8//Crf, the steady state corresponding to x* 
will be locally stable. Once condition ( l23l is violated, stability 
of x* is lost via a supercritical Hopf bifurcation and a limit cycle 
is born. One may even compute the Hopf period of this limit 
cycle by assuming that A = ja>H (j = V— T) in Equation l22l 
where <dh is the Hopf angular frequency. Equating real and 
imaginary parts of the resultant yields u>h = y/a^/ai or 



2tt 



T„ = — = 2n x 



d-KdfOnuyi 



These local stability results tell us nothing about the global 
behaviour when stability is lost, but it is possible to characterize 
the global behaviour of a repressible operon with the following 



Theorem 2. ASmitn 179951 Theorem 4.1 & Theorem 4.2, Chap- 
ter 3) For a repressible operon with ip given by Equation^ de- 
fine I r = [Ki/K, 1]. There is a globally attracting box Br C 
defined by 

B R = {(x u x2,x 3 ) : Xi 6 I R , i= 1,2,3} 



such that the flow S is directed inward everywhere on the sur- 
face of Br. Furthermore there is a single steady state X* e Br. 
If X* is locally stable it is globally stable, but if X* is unsta- 
ble then a ge neralization of the Poincare-Bendixson theorem 
\Smith\\199?\ Chapter 3) implies the existence of a globally sta- 
ble limit cycle in Br. 



Remark 2. There is no necessary connection between the Hopf 
period computed from the local stability analysis and the period 
of the globally stable limit cycle. 

3. Fast and slow variables 

In dynamical systems, considerable simplification and in- 
sight into the behaviour can be obtained by identifying fast and 
slow variables. This technique is especially useful when one 
is initially interested in the approach to a steady state. In this 
context a fast variable is one that rela xes much mor e rapidly to 
an equilibrium than a slow variable (IHakenl 119831) . In many 
systems, including chemical and biochemical ones, this is of- 
ten a consequence of differences in degradation rates, with the 



6 



fastest variable the one that has the largest degradation rate. We 
employ the same strategy here to obtain approximations to the 
population level dynamics that will be used in the next section. 

It is often the case that the degradation rate of mRNA is much 
greater than the corresponding degradation rates for both the 
intermediate protein and the effector (71 » 72, 73) so in this 
case the mRNA dynamics are fast and we have the approximate 
relationship 

Xi a K d f{x-$). 

Consequently the three variable system describing the generic 
operon reduces to a two variable one involving the slower inter- 
mediate and effector: 



-IT = 72[*rf/(*3) - X2], 
at 

dx 3 

—r- = 73(*2 - X 3 ). 

at 



(24) 
(25) 



In our considerations of specific single operon dynamics be- 
low we will also have occasion to examine two further subcases, 
namely 

Case 1. Intermediate (protein) dominated dynamics. If it 

should happen that j\ » 73 » 72 (as for the lac operon, then 
the effector also qualifies as a fast variable so 

X3 - x 2 

and thus from (I24b -d25l> we recover the one dimensional equa- 
tion for the slowest variable, the intermediate: 



dx2 

—77 = 72[Kdf(X2) - x 2 ]. 
at 



(26) 



Case 2. Effector (enzyme) dominated dynamics. Alternately, 
if 71 » 72 » 73 then the intermediate is a fast variable relative 
to the effector and we have 



x 2 x 3 

so our two variable system ( |241-(|25]> reduces to a one dimen- 
sional system 

-77 = 73 [Kdf(x3) - x 3 ] (27) 
at 

for the relatively slow effector dynamics. 
Both Equations l26l and 1271 are of the form 



dx 
dt 



y[Kdf{x) - x] 



(28) 



where 7 is either 72 for protein (x 2 ) dominated dynamics or 73 
for effector (^3) dominated dynamics. 

4. Distributions with intrinsic bursting 

4.1. Generalities 



It is wel l documented experim e ntally ( Cai et 



Sigal et all 120061: lYu et al 



Chubb et all l2006t iGolding et all 120051: iRai et 
: all 120061 \\ 



( Cai etall 
IRai et all 



2006 



2006 



20061) that in some organisms the 



amplitude of protein production through bursting translation of 



mRNA is exponentially distributed at the single cell level with 
density 

h(y) = le- yrb , (29) 



where b is the average burst size, and that the frequency of 
bursting ip is dependent on the level of the effector. Writing 
Equation[29]in terms of our dimensionless variables we have 



h(x) = -e- x/b . 
b 



(30) 



Remark 3. The technique of eliminating fast variables de- 
scribed in Section \23\above (als o known as the adiabatic elim- 
ination technique Uiakem \1983\) ) has been extended to stochas- 
tically perturbed systems whe n the perturbation is a Gaussian 



distributed wh ite noise, c.f. \ Stratonovich\ (179631 Chapter 4, 



Sectio n 11.1), Wilemski 1 197a) . Titular i l978\) . and Gardiner 



U983l Section 6.4). However, to the best of our knowledge, this 
type of approximation has never been extended to the situation 
dealt with here in which the perturbation is a jump Markov pro- 
cess. 

The single cell analog of the population level intermediate 
protein dominated Case 1 above (when 71 » 73 » 72) is 



dxi 
dt 



= -Jixi + a(h, tf(x 2 )), with <p(x 2 ) = y2iPmf{x 2 ), (31) 



where 3(h, <p) denotes a jump Markov process, occurring at a 
rate ip, whose amplitude is distributed with density h as given in 
( 1301 . Analogously, in the Case 2 effector dominated situation 
the single cell equation becomes 



dxi 
dt 



= -73x3 + H(/2, tp(x 3 )), with ip(x 3 ) = jiip m f{x 3 ). (32) 



Equations [3T1andl32lcan both be written as 

dx 

— = -y X + Z(h,ip(x)), with p(x) = yKbf(x), K b = <p m . 
dt 

Remark 4. In the case of bursting we will always take Kb = tp„ 
in contrast to the deterministic case where kj — b,np m . 



From iMackev and Tvran-Kamiriskal d2008l) the correspond- 
ing operator equation for the evolution of the density u(t,x) 
when there is a single dominant slow variable is given by 

du(t,x) d(xu(t,x)) 

- 7 = -yK b f(x)u(t, x) 



dt 



dx 



7 Kb \ 

Jo 



f(y)u(t,y)h(x-y)dy. 



(33) 



Re mark 5 . This is a straightforward generalization of what 
\Gardiner\ \l983 Section 3.4) refers to as the differential 
Chapman-Kolmogorov equation. 

Stationary solutions u*(x) of d3~3l are solutions of 

d(xu t (x)) 



dx 



■I 



K b f(x)u«{x)+K b f(y)u t (y)h(x-y)dy. (34) 



If there is a unique stationary density, then the solution 
u(t,x) of Equation 1331 is said to be asymptotically stable 
W9 



dLasota and Mackevill994l) in the sense that 

I \u(t, x) - u t (x)\dx = 
Jo 



lim 

t— >oo 



for all initial densities u(0, x). 



Theorem 3. wlackev and T/ran-K amihska 2008 Theorem 7). 
The unique stationary density of Equation\34\ with f given by 
Equation\9\and h given by ( I29H . is 



, , c 

u t (x) = —e 



-x/b 



exp 



K b dy 

J y 



(35) 



where C is a normalizing constant such that J Q u*(x)dx — 1. 
Further, u(t, x) is asymptotically stable. 

Remark 6. The stationary density ( 1551 ) is found by rewriting 
Equation\34\in the form 

^j X ^ + ^y-^ - Kt^^-y(x) = 0, y(x) = xu„(x) 
dx b x 

using Laplace transforms and solving by quadratures. Note 
also that we can represent u* as 



u t (x) = C exp 



Kbfiy) l n 

y by)' 



where C is a normalizing constant. 

4.2. Distributions in the presence of bursting 
4.2.1. Protein distribution in the absence of control 

If the burst frequency ip - yn\,f is independent of the level 
of all of the participating molecular species, then the solution 
given in Equation[35]is the density of the gamma distribution: 



w»(x) = 



1 



b K "T(K b ) 



where T(-) denotes the gamma function. For k/, e (0, 1), m»(0) = 
oo and u, is decreasing while for Kb > 1, u„(Q) = and there is 
a maximum at x = b(K\, - 1). 

4.2.2. Controlled bursting 

We next consider the situation in which the burst frequency 
<p is dependent on the level of x, c.f. Equation|5] This requires 
that we evaluate 

where A, A are enumerated in Table[T]for both the inducible and 
repressible operons treated in Section l272l and 



nA \ A 



Consequently, the steady state density d35l l explicitly becomes 




Figure 5: Schematic illustration of the possibility of one, two or three solutions 
of Equation f38]f or varying values of ki, with bursting inducible regulation. The 
straight lines correspond (in a clockwise direction) to ki, e (Q,Kb-), ki, = Kb-, 
ki, e (ati,_, Kf,+) (and respectively k\, < K, Kb = K, K < k/,), Kb = Kb+, and 
Kb+ < Kb- This figure was constructed with n = 4, K = 10 and b = 1 for which 
ki,- = 4.29 and ki,+ = 14.35 as computed from )42t . See the text for further 
details. 



The first two terms of Equation[36]are simply proportional to 
the density of the gamma distribution. For < Kb AT 1 < 1 we 
have m,(0) = oo while for /q,A 1 > 1, «*(0) = and there is at 
least one maximum at a value of x > 0. We have u*{x) > for 
all x > and from Remarket follows that 



u'„(x) = u*(x) 



K b f(x) 



x>Q. 



(37) 



u»{x) = Ce 



~x/b« b A- ] -\ 



(A + Ax") 6 



(36) 



Observe that if Kb < 1 then u» is a monotone decreasing function 
of x, since Kbf(x) < 1 for all x > 0. Thus we assume in what 
follows that ki, > 1 . 

Since the analysis of the qualitative nature of the stationary 
density leads to different conclusions for the inducible and re- 
pressible operon cases, we consider each in turn. 

4.2.3. Bursting in the inducible operon 

For 6 > 0, as in the case of an inducible operon, the third 
term of Equation[36]is a monotone increasing function of x and, 
consequently, there is the possibility that u t may have more than 
one maximum, indicative of the existence of bistable behaviour. 
In this case, the stationary density becomes 

«.(*) = Ce-*l b x*» K ~ l -\K + x n ) e , 8 = % - K- 1 ). 

n 

From (l37l i it follows that we have u' t (x) = for x > if and 
only if 

1 Ix \ 1 + x" 

— - + 1 = . (38) 

K b \b ) K + x" 

Again, graphical arguments (see Figure |5]l show that there may 
be up to three roots of d38l . For illustrative values of n, K, and 
b, Figure|6]shows the graph of the values of x at which u' t {x) = 



8 



as a function of Kb- When there are three roots of (1381 1. we label 
them as x\ < X2 < A3. 

Generally we cannot determine when there are three roots. 
However, we can determine when there are only two roots x\ < 
A3 from the argument of Section l2.3.2l At x\ and A3 we will not 
only have Equation[38] satisfied but the graph of the right hand 
side of (l38T l will be tangent to the graph of the left hand side at 
one of them so the slopes will be equal. Differentiation of ( l38l l 
yields the second condition 



A'" 



1 



(K + x") 2 K h b(K - 1) 



(39) 



We first show that there is an open set of parameters (b, K, Kb) 
for which the stationary density u, is bimodal. From Equations 
[38]and[39]it follows that the value of x+ at which tangency will 
occur is given by 

x± = b(K b - l)z± 
and z± are positive solutions of equation 

£ = l_ z _^l_ z) 2, where B = 

n (K-V)K b 



We explicitly have 



z ± = ^-\2Bn-{n + \). 
2Bn ' 



+ i f - ABn 



provided that 



(n+1) 2 >/3 K(K h -l) 



(40) 



An (K-l)K b ' 
Equation l40l is always satisfied when k\, < K or when Kb > K 
and K is as in the deterministic case ( fT9b . Observe also that we 
have z+ > > Z- for Kb < K and z+ > Z- > for Kb > K. The 
two corresponding values of b at which a tangency occurs are 
given by 



b ± = 



1 



K 



(K b -l)z ± V/3(l -z ± ) 



- K and z+ > 0. 



If ki, < K then m»(0) = 00 and w» is decreasing for b < b+, 
while for b > b+ there is a local maximum at x > 0. If Kb > K 
then m*(0) = and u* has one or two local maximum. As a 
consequence, for n > 1 we have a bimodal steady state density 
u* if and only if the parameters Kb and K satisfy f40) . Kb > K, 
and b e (b+,bJ). 

We now want to find the analogy between the bistable be- 
havior in the deterministic system and the existence of bimodal 
stationary density w,. To this end we fix the parameters b > 
and K > 1 and vary Kb as in Figure|5] Equations I38landl39lcan 
also be combined to give an implicit equation for the value of 
x± at which tangency will occur 



x 2 " -(K-l) 



K+l 



K- 1 



x" - nb(K ■ 



l)x"- 1 + K ■■ 



and the corresponding values of ki,± are given by 
<x T + b\(K + xl' 



«b± = 



1+xZ 



(41) 



(42) 



There are two cases to distinguish. 

Case 1. < Kb < K. In this case, m„(0) = 00. Further, the same 
graphical considerations as in the deterministic case show that 
there can be none, one, or two positive solutions to Equation 
[381 If ki, < Kb-, there are no positive solutions, u* is a monotone 
decreasing function of x. If Kb > Kb-, there are two positive 
solutions (x2 and I3 in our previous notation, x\ has become 
negative and not of importance) and there will be a maximum 
in m* at X3 with a minimum in u, at X2- 

Case 2. < K < Kb. Now, m»(0) = and there may be one, 
two, or three positive roots of Equation[38] We are interested in 
knowing when there are three which we label as x\ < iS < %3 
as Xi,X3 will correspond to the location of maxima in u„ while 
X2 will be the location of the minimum between them and the 
condition for the existence of three roots is k\,- < Kb < ki,+. 

We see then that the different possibilities depend on the re- 
spective values of K, Kb-, Kb+, and Kb- To summarize, we may 
characterize the stationary density m* for an inducible operon in 
the following way: 

1. Unimodal type 1: u*(0) = 00 and u» is decreasing for 
< Kb < Kb- and < Kb < K 

2. Unimodal type 2: m»(0) = and u* has a single maxi- 
mum at 

(a) x\ > for K < Kb < Kb- or 

(b) at A3 > for k/,+ < Kb and K < Kb 

3. Bimodal type 1: m»(0) = 00 and u* has a single maximum 
at A3 > for Kb- < Kb < K 

4. Bimodal type 2: m*(0) = and m* has two maxima at 
Ai, A3, < x\ < A3 for Kb- < ki, < Kb+ and K < Kb 

Remark 7. Remember that the case n — 1 cannot display testa- 
bility in the deterministic case. However, in the case of bursting 
in the inducible system when n — 1, if j + 1 < Kb < K and 
b > -^j, then m„(0) = 00 and u* also has a maximum at A3 > 0. 
Thus in this case one can have a Bimodal type 1 stationary den- 
sity. 

We now choose to see how the average burst size b affects 
bistability in the density u* by looking at the parametric plot of 
Kb(x) versus K(x). Define 



F(x,b) = 



a" + 1 



nx"- l (x + b) 



(43) 



Then 



I + y p( x h\ x + 

K(x,b) = - and K b (x,b) = [K(x,b)+x n ]- 



1 - F(x, b) 



b(x" + 1) 
(44) 

The bifurcation diagram obtained from a parametric plot of K 
versus Kb (with a as the parameter) is illustrated in Figure [7] 
for n — 4 and two values of b. Note that it is necessary for 
< K < Kj, in order to obtain Bimodal type 2 behaviour. 

For bursting behaviour in an inducible situation, there are two 
different bifurcation patterns that are possible. The two differ- 
ent cases are delineated by the respective values of K and Kb, 
as shown in Figure [6] and Figure [7] Both bifurcation scenarios 
share the property that while increasing the bifurcation param- 
eter Kb from to 00, the stationary density m, passes from a 



9 




Figure 6: Full logarithmic plot of the values of x at which u' t (x) = versus the 
parameter Kb. obtained from Equation [38] for n = 4, K = 10, and (left to right) 
b = 5, 1 and b = . Though somewhat obscured by the logarithmic scale for x, 
the graphs always intersect the Kb axis at fq, = K. Additionally, it is important 
to note that «»(0) = for K < Kb, and that there is always a maximum at for 
< Kb < K. See the text for further details. 




5 10 15 

K 

Figure 7: In this figure we present two bifurcation diagrams (for n = 4) in 
(K, Kt>) parameter space delineating unimodal from bimodal stationary densities 
«» in an inducible operon with bursting as obtained from Equations l44l with l43l 
The upper cone-shaped plot is for b = while the bottom one is for b = 1 . In 
both cone shaped regions, for any situation in which the lower branch is above 
the line Kb = K (lower straight line) then bimodal behaviour in the stationary 
solution h*( y) will be observed with maxima in u, at positive values of x, x\ 
and x-j. 



unimodal density with a peak at a low value (either or xi) to 
a bimodal density and then back to a unimodal density with a 
peak at a high value (X3). 

In what will be referred as Bifurcation type 1, the maximum 
at x — disappears when there is a second peak at x = X3. The 
sequence of densities encountered for increasing values of Kb is 
then: Unimodal type 1 to a Bimodal type 1 to a Bimodal type 2 
and finally to a Unimodal type 2 density. 

In the Bifurcation type 2 situation, the sequence of density 
types for increasing values of Kb is: Unimodal type 1 to a Uni- 
modal type 2 and then a Bimodal type 2 ending in a Unimodal 
type 2 density. 

The two different kinds of bifurcation that can occur are eas- 
ily illustrated for b — 1 as the parameter Kb is increased. (An 
enlarged diagram in the region of interest is shown in Figure[8]) 
Figure|9]illustrates Bifurcation type 1, when K — 4, and Kb in- 
creases from low to high values. As Kb increases, we pass from 
a Unimodal type 1 density, to a Bimodal type 1 density. Further 
increases in Kb lead to a Bimodal type 2 density and finally to 
a Unimodal type 2 density. This bifurcation cannot occur, for 
example, when b — 4j and K < 15 (see Figure [7) . 

Figure [10] shows Bifurcation type 2, when K = 3. As Kb in- 
creases, we pass from a Unimodal type 1 density, to a Unimodal 
type 2 density. Then with further increases in Kb, we pass to a 
Bimodal type 2 density and finally back to a Unimodal type 2 
density. 

Remark 8. There are several qualitative conclusions to be 
drawn from the analysis of this section. 

1 . The presence of bursting can drastically alter the regions 
of parameter space in which bistability can occur relative 
to the deterministic case. Figure \TT\presents the regions 
of bistability in the presence of bursting in the (K,b ■ Kb) 
parameter space, which should be compared to the region 
of bistability in the deterministic case in the (K, /q/) param- 
eter space (bKb is the mean number of proteins produced 
per unit of time, as is Kd) 
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Figure 9: In this figure we illustrate Bifurcation type 1 when intrinsic bursting 
is present. For a variety of values of the bifurcation parameter k\, (between 3 
and 6 from top to down), the stationary density »* is plotted versus x between 
and 8. The values of the parameters used in this figure are b = 1, K = 4, and 
n = 4. For k/, < 3.5, u„ has a single maximum at x = 0. For 3.5 < Kj, < 4, 
n» has two local maxima at x = and I3 > 1. For 4 < Kb < 5.9, n» has two 
local maxima at < .?i < A3. Finally, for /q, > 5.9, u„ has a single maximum 
at X3 > 1. Note that for each plot of the density, the scale of the ordinate is 
arbitrary to improve the visualization. 



2. When < Kb < K, at a fixed value of Kb, increasing the av- 
erage burst size b can lead to a bifurcation from Unimodal 
type 1 to Bimodal type 1. 

3. When < K < Kb, at a fixed value of Kb, increasing b can 
lead to a bifurcation from Unimodal type 2 to Bimodal type 
2 and then back to Unimodal type 2. 

4.2.4. Bursting in the repressible operon 

The possible behaviours in the stationary density u t for the 
repressible operon are easy to delineate based on the analysis 
of the previous section, with Equation[38]replaced by 



_j_/x \ _ 1 +x" 
K b \b + / 1 + Ax" ' 



(45) 



Again graphical arguments (see Figure ITZb show that Equation 
|45]may have either none or one solution. Namely, 

1. For < Kb < 1, m*(0) = 00 and u, is decreasing. Equation 
|45]does not have any solution (Unimodal type 1). 

2. For 1 < Kb, m*(0) = and u* has a single maximum at a 
value of x > determined by the single positive solution 
of Equationl45l(Unimodal type 2). 

4.3. Recovering the deterministic case 

We can recover the deterministic behaviour from the bursting 
dynamics with a suitable scaling of the parameters and limiting 
procedure. With bursting production there are two important 
parameters (the frequency k/, and the amplitude b), while with 
deterministic production there is only kj. The natural limit to 
consider is when 

b — > 0, /<•/,—> 00 with bKb = Kd- 



Figure 10: An illustration of Bifurcation type 2 for intrinsic bursting. For 
several values of the bifurcation parameter k\, (between 2.8 and 5 from top 
to down), the stationary density u„ is plotted versus x between and 8. The 
parameters used are b = 1, K = 3, and n = 4. For ki, < 3, u, has a single 
maximum at x = 0, and for 3 < k/, < 3.3, u, has a single maximum at .?i > 0. 
For 3.3 < Kj, < 4.45, w* has two local maxima at < .?i < X3, and finally for 
K b ^ 4.45 u„ has a single maximum at X3 > 0. Note that for each plot of the 
density, the scale of the ordinate is abritrary to improve the visualization. 



In this limit, the implicit equations which define the maximum 
points of the steady state density, become the implicit equations 
( TT3T > and ( fTo*b which define the stable steady states in the deter- 
ministic case. 

The bifurcations will also take place at the same points, be- 
cause we recover Equation [18] in the limit. However, Bimodal- 
ity type 1 as well as the Unimodal type 1 behaviours will no 
longer be present, as in the deterministic case, because for 
ki, > co we have k\, > K. Finally, from the analytical ex- 
pression for the steady-state density (f36t m* will became more 
sharply peaked as b — > 0. Due to the normalization constant 
(which depends on b and Kb), the mass will be more concen- 
trated around the larger maximum of u*. 

5. Distributions with fluctuations in the degradation rate 

5.1. Generalities 

For a generic one dimensional stochastic differential equation 
of the form 

dx(t) — a(x)dt + o~(x)dw(t) 
the corresponding Fokker Planck equation 



du d(au) 1 3 2 (ct 2 m) 
dt dx 2 dx 2 

can be written in the form of a conservation equation 

du ^ dJ 
dt dx 



(46) 



where 



1 d(o- 2 u) 
J — a u 



2 dx 

is the probability current. In a steady state when d,u = 0, the 
current must satisfy / = constant throughout the domain of 
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Figure 1 1 : The presence of bursting can drastically alter regions of bimodal 
behaviour as shown in this parametric plot (for n = 4) of the boundary in (K, b ■ 
Kl,) parameter space delineating unimodal from bimodal stationary densities u„ 
in an inducible operon with bursting and in (K.icj) parameter space delineating 
one from three steady states in the deterministic inducible operon. From top to 
bottom, the regions are for b = 10, b = 1, b = 0.1 and b = 0.01. The lowest 
(heavy dashed line) is for the deterministic case. Note that for b — 0.1, the two 
regions of bistability and bimodality coincide and are indistinguishable from 
one another. 



the problem. In the particular case when J — at one of the 
boundaries (a reflecting boundary) then 7 = for all x in the 
domain and the steady state solution u* of Equation|46]is easily 
obtained with a single quadrature as 



m*(x) = 



^) 6XP { 2 f 



a(y) 



dy\, 



0~ 2 (x) \ J cr 2 (y) 
where C is a normalizing constant as before. 



5.2. Fluctuations in degradation rate 

In our considerations of the effects of continuous fluctua- 
tions, we examine the situation in which fluctuations appear 
in the degradation rate y of the gen eric equation (l28l . From 



standard chemical kinetic arguments (lOppenheim et al. , 1969) 



if the fluctuations are Gaussian distributed the mean numbers of 
molecules decaying in a time dt is simply yxdt and the standard 
deviation of these numbers is proportional to -\fx. Thus we take 
the decay to be given by the sum of a deterministic component 
yxdt and a stochastic component cr -\fxdw(t), where w is a stan- 
dard Brownian motion, and write Equation [28] as a stochastic 
differential equation in the form 

dx = y[i<df(x) - x]dt + cr yfxdw. 

Within the Ito interpretation of stochastic integration, this 
equation has a corresponding Fokker Planck equation for 
the evolution of the en semble density u(t,x) given by 
dLasota and Mackevi[l994l) 



du d [(yi<df(x) — yx)u] cr 2 d (xu) 
dt dx 2 dx 2 



(47) 



Figure 12: Schematic illustration that there can be one or no solution of Equa- 
tion !45l depending on the value of aj, with repressible regulation. The straight 
lines correspond (in a clockwise direction) to k\, = 2 and k\, = 0.8. This figure 
was constructed with n = 4, A = 10 and b = 1. See the text for further details. 



In the situation we consider here, cr(x) = cr y[x and a(x) = 
yKdf(x) - yx. Further, since concentrations of molecules cannot 
become negative the boundary at x — is reflecting and the 
stationary solution of Equationl47lis given by 



< ^ C 
u*(x) = —e 



-2yx/o~ 2 



exp 



2jK d 



ff 



dy 



Set K e = lyKdlcr 2 . Then the steady state solution is given ex- 
plicitly by 

u,(x) = Ce- 3 * x,ai 3«*~ l - 1 [A + Ax"] 6 , (48) 
where A, A > and 9 are given in Table Q] 
Remark 9. Two comments are in order. 

1 . Because the form of the solutions for the situation with 
bursting (intrinsic noise) and extrinsic noise are identical, 
all of the results of the previous section can be carried 
over here with the proviso that one replaces the average 
burst amplitude b with b — ■> cr 2 /2y = b K and /(/,—> K e — 
2yK d lcr 2 = Kd/K- 

2. We can look for the regions of bimodality in the (K, k c i)- 
plane, for a fixed value ofb w . We have the implicit equa- 
tion for x ± 



(K-l) 



K+ 1 



K- 1 



x" - nb w (K - l)x"- 1 +K = 



and the corresponding values of kj are given by 

'K + x'l' 



Kd± = (x* + b w ) 



1 + x% 



Then the bimodality region in the (K, Kd)-plane with noise 
in the degradation rate is the same as the bimodality region 
for bursting in the (K, bKt)-plane. 
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We have also the following result. 
Theorem 4. ( Pichor and Rudnicki 200(\ Theorem 2). The 



unique stationary density of Equation^7\is given by Equation 
\48\ Further u(t, x) is asymptotically stable. 

5.3. The deterministic limit 

Here again we can recover the deterministic behavior from a 
limit in the extrinsic fluctuations dynamics. In this case, how- 
ever, the frequency and the amplitude of the perturbation are 
already scaled. Then the limit cr — » gives the same result as 
in the deterministic case. 



6. Discussion and conclusions 

In trying to understand experimentally observed distribu- 
tions of intracellular components from a modeling perspec- 
tive, the norm in computational and syst ems biology is of ten 
to use algorithms developed initially by iGillespid d 19771) to 
solve the chemical ma ster equation for specific situations. See 
Lipniacki et al. (2006) for a typical example. However these in- 



vestigations demand long computer runs, are computationally 
expensive, and further offer little insight into the possible diver- 
sity of behaviours that different gene regulatory networks are 
capable of. 

There have been notable exceptions in which the prob- 
lem ha s been treated fr o m an an alytical point of view , 



c.f. Kepler and Elston (2001), Friedman et 



Bobrowski et al 



"aD 
uinl 



(2006), 
d2008al) . 



(120071) . and IShahreza ei and Swai 
The advantage of an analytic development is that one can de- 
termine how different elements of the dynamics shape temporal 
and steady state results for the densities u(t, x) and u,{x) respec- 
tively. 

Here we have extended this analytic treatment to simple sit- 
uations in which there is either bursting transcription and/or 
translation (building o n and expanding the original work of 
( Friedman et al. , 2006|)). or fluctuatio ns in degradation rates, as 



an alternative to the IGillespid d 19771) algorithm approach. The 
advantage of the analytic approach that we have taken is that it 
is possible, in some circumstances, to give precise conditions on 
the statistical stability of various dynamics. Even when analytic 
solutions are not available for the partial integro-differential 
equations governing the density evolution, the numerical solu- 
tion of these e quations may be computationally more tractable 
than using the IGillespid d 1977b approach. 

One of the more surprising results of the work reported here 
is that the stationary densities in the presence of bursting noise 
derived in Section [4] are analytically indistinguishable from 
those in the presence of degradation noise studied in Section 
|5] We had expected that there would be clear differences that 
would offer some guidance for the interpretation of experimen- 
tal data to determine whether one or the other was of predomi- 
nant importance. Of course, the next obvious step is to examine 
the problem in the presence of both noise sources simultane- 
ously. However, the derivatio n of the evolution equation in this 
case, as has been pointed out dHierro and Dopazo , l2009h . is not 



straightforward and we will report on our results in a separate 
communication. 

In terms of the issue of when bistability, or a unimodal versus 
bimodal stationary density is to be expected, we have pointed 
out the analogy between the unimodal and bistable behaviours 
in the deterministic system and the existence of bimodal sta- 
tionary densities in the stochastic systems. Our analysis makes 
clear the critical role of the dimensionless parameters n, k (be 
it Kd, Kb, or K e ), b (either b or b w ), and the fractional leakage 
K . The relations between these defining the various possible 
behaviours are subtle, and we have given these in the relevant 
sections of our analysis. 

The appearance of both unimodal and bimodal distri- 
butions of molecular constituents as well as what we 
have termed Bifurcation Type 1 and Bifurcation Type 
2 have been extensively discussed in the applie d math - 



ematics literature (c.f. Horsthe mke and Lefeveri d!984l) . 



Feistel and Ebelinp dl989h and others) and the bare founda- 



tio ns of a s tochas tic bifurcation theory have been laid down 
by I Arnold! J 19981) . Significantly, these are also well doc- 



bv Gardner et al. (2000b. Acaretal.1 (2005b. Friedman et al. 


(2006 


). Hawkins and Smolke ( 20061). IZac 


larioudakis et al. 


(2007 


). Mariani et all (20101). and ISong et al. 


d2010h for both 



prokaryotes and eukaryotes. If the biochemical details of a par- 
ticular system are sufficiently well characterized from a quan- 
titative point of view so that relevant parameters can be es- 
timated, it may be possible to discriminate between whether 
these behaviours are due to the presence of bursting transcrip- 
tion/translation or extrinsic noise. 
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